#include "TString.h"
#include "TH2D.h"
#include "TTree.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLine.h"

#include "TLegend.h"
#include "style.h"
#include "NeutrinoTools.h"


/*
 enum eventtype{
    kZero = 0,
    kPionPlus=1,
    kPionMinus=2,
    kProton = 3,
    kOther = 4,
    kHasParent = 5,

    kMPion = 6,
    kMProton = 7,
    k0Muon = 8,
    k0Proton = 9,
    kMMuon = 10,
    kBad = 11
  };

 enum pid{
    kMuon=0,
    kProton=1,
    kPion=2,
    kElectron=3,
    kOther=4
  };

*/


TString CCQEcut()
{
  return "fSelX_Charge==-999";
}

TString CCQEPIDcut()
{
  return "&& fProtonSimPDG == 2212 && fMuonSimPDG == 13";
}

TString CCQETruePDGModeAlltypeNoECalcut(TString & alltit)
{
  alltit.Prepend("true PID, CCQE, N^{track}_{ECal}=0");
  return CCQEPIDcut()+" && fNeutMode==1 && fNECalTracks==0";
}

TString CCQETruePDGModeAlltypecut(TString & alltit)
{
  alltit.Prepend("true PID, CCQE");
  return CCQEPIDcut()+" && fNeutMode==1";
}

TString CCQETruePDGModeTypecut(TString & alltit)
{
  alltit.Prepend("true PID, CCQE, no nuclear emission");
  return CCQEPIDcut()+" && fNeutMode==1 && fEventTypeSim==0";
}

TString CCQETruePDGModeBreakUpcut(TString & alltit)
{
  alltit.Prepend("true PID, CCQE, w/ nuclear emission");
  return CCQEPIDcut()+" && fNeutMode==1 && fEventTypeSim!=0";
}

TString CCPXcut()
{
  return "fSelX_Charge==1 && fSelX_PID == 2";
}

TString CCMuPPcut()
{
  return "fSelX_Charge==1 && fSelX_PID == 1";
}

TString CCPXWithECalcut()
{
  return "fSelX_Charge==1 && fSelX_PID == 2 && fNECalTracks==0";
}

TString CCPXWrongChargeBgcut(TString & alltit)
{
  alltit.Prepend("#pi^{-}");
  return "fSelX_Charge==-1 && fSelX_PID == 2";
}

TString DeltaPIDcut()
{
  return "&& fXSimPDG == 211 && fProtonSimPDG == 2212 && fMuonSimPDG == 13";
}

TString MuPPPIDcut()
{
  return "&& fXSimPDG == 2212 && fProtonSimPDG == 2212 && fMuonSimPDG == 13";
}

TString CCMuPPTruePDGTypecut(TString & alltit)
{
  alltit.Prepend("true PID, #mu pp, no nuclear emission");
  return MuPPPIDcut()+" && fEventTypeSim==3";
}

TString CCPXTruePDGModeTypecut(TString & alltit)
{
  alltit.Prepend("true PID, #Delta^{++}, no nuclear emission");
  return DeltaPIDcut()+" && fNeutMode==11 && fEventTypeSim==1";
}

TString CCPXTruePDGModeBreakUpcut(TString & alltit)
{
  alltit.Prepend("true PID, #Delta^{++}, w/ nuclear emission");
  //&& fTopology!=1 very small Neut 11, Genie 66
  return DeltaPIDcut()+" && fNeutMode==11 && fEventTypeSim!=1";
}

TString CCPXTruePDGModeAlltypecut(TString & alltit)
{
  alltit.Prepend("true PID, #Delta^{++}");
  return DeltaPIDcut()+" && fNeutMode==11";
}

TString CCPXTruePDGBgcut(TString & alltit)
{
  alltit.Prepend("true PID, non-#Delta^{++} bg");
  return DeltaPIDcut()+" && fNeutMode!=11";
}

void SetTagTitle(const TString fin, TString &tag, TString &alltit)
{
  if(fin.Contains("Neut")){
    tag+="Neut";
    alltit.Prepend("Neut ");
  }
  else if(fin.Contains("Genie")){
    tag+="Genie";
    alltit.Prepend("Genie ");
  }

  if(fin.Contains("WaterOnly")){
    tag+="Water";
    alltit.Prepend("Water ");
  }
  else if(fin.Contains("AirOnly")){
    tag+="Air";
    alltit.Prepend("Air ");
  }
}

void plotDoubleT(const TString fin, const Int_t opt=1)
{
  TFile::Open(fin);

  style::SetGlobalStyle();

  TTree * tree= (TTree*) gDirectory->Get("tree");

  TString cut, chan, tag;

  if(opt==1){
    cut=CCQEcut();
    chan = " && fEventTypeSim==0"; 
    tag="exclusive #mu^{-}p event";
  }
  else if(opt==2){
    cut=CCQEcut();
    chan = " && fNeutMode==1"; 
    tag="CCQE";
  }
  else if(opt==3){
    cut=CCQEcut();
    chan = " && fTopology==0"; 
    tag="CC0pi";
  }
  else if(opt==4){
    cut=CCPXcut();
    chan = " && fNeutMode==11"; 
    tag="#Delta^{++}";
  }
  else if(opt==5){
    cut=CCPXcut();
    chan = " && fEventTypeSim==1"; 
    tag="exclusive #mu^{-}p#pi^{+} event";
  }
  else if(opt==6){
    cut=CCPXcut();
    chan = " && fTopology==1"; 
    tag="CC1pi";
  }
  else if(opt==7){
    cut=CCPXcut();
    chan = " && fNeutMode==11 && fEventTypeSim==1"; 
    tag="exclusive #mu^{-}p#pi^{+} event from #Delta^{++}";
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  TCanvas * c1 =new TCanvas;
  style::PadSetup(c1);

  TH1D * h0=new TH1D("h0","",5,0,5);
  TH1D * h1=new TH1D("h1","",5,0,5);

  h1->SetTitle(Form(";N_{ECal}^{track}; %s purity", tag.Data()));

  tree->Draw("fNECalTracks>>h0", cut);
  tree->Draw("fNECalTracks>>h1", cut + chan,"same");
  TF1 * ff= new TF1("ff","[0]",0,10); 
  ff->SetParameter(0,h1->GetEntries()/h0->GetEntries()); 
  h1->Divide(h0); 
  style::ResetStyle(h1);
  h1->SetMinimum(0);
  h1->Draw(); 
  style::ResetStyle(ff);
  ff->Draw("same");

  tag.ReplaceAll("#","");
  tag.ReplaceAll("+","");
  tag.ReplaceAll("{","");
  tag.ReplaceAll("}","");
  tag.ReplaceAll("^","");
  tag.ReplaceAll(" ","");

  c1->Print(Form("outplot/%spurity.png",tag.Data()));
}

void plotVertex(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();

  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);

  TString cut, tag, topo;
  if(opt==0){
    cut=CCQEcut();
    topo="&& fTopology==0";
    tag="cc0pi";
  }
  else if(opt==1){
    cut=CCPXcut();
    topo="&& fTopology==1";
    tag="cc1pi";
  }
  else if(opt==2){
    cut=CCPXcut();
    topo="&& fNeutMode==11";
    tag="delta";
  }
  else if(opt==3){
    cut=CCQEcut();
    topo="&& fNeutMode==1";
    tag="ccqe";
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  TString alltit;
  SetTagTitle(fin, tag, alltit);
  
  tt->Draw("fTopology", cut);
  cc->Print(Form("outplot/%sfTopology.png", tag.Data()));  

  tt->Draw("fNeutMode", cut);
  cc->Print(Form("outplot/%sfNeutMode.png", tag.Data()));

  tt->Draw("fEventTypeSim", cut);
  cc->Print(Form("outplot/%sfEventTypeSim.png", tag.Data()));

  gPad->SetLogy();

  //2nd proton and multi proton
  tt->Draw("fEventTypeSim", cut+topo);
  cc->Print(Form("outplot/%stopo_fEventTypeSim.png", tag.Data()));

  //gamma and neutron
  tt->Draw("fSimOtherPDG", cut+topo+" && fEventTypeSim==4 && (fSimOtherPDG!=22)");
  cc->Print(Form("outplot/%stopo_other_nogamma_pdg.png", tag.Data()));



  //gamma and neutron and
  Double_t ent= tt->Draw("fSimOtherPDG", cut+topo+" && fEventTypeSim==4 && (fSimOtherPDG!=22 && fSimOtherPDG!=2112)");

  if(ent==0){
    return;
  }

  cc->Print(Form("outplot/%stopo_other_nogamma_noneutron_pdg.png", tag.Data()));
}

void plotECalPDG(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();

  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);

  TString cut, tag, topo, var, truecut;
  if(opt%10==0){
    cut=CCQEcut()+" && fProtonSimPDG == 2212 && fMuonSimPDG == 13";
    tag="ccqe";
  }
  else if(opt%10==1){
    cut=CCPXcut()+" && fProtonSimPDG == 2212 && fMuonSimPDG == 13";
    tag="ccpx";
  }
  else if(opt%10==2){
    cut=CCQEcut()+" && fProtonSimPDG == 2212 && fMuonSimPDG == 13 && fNeutMode==1 && fEventTypeSim ==0";
    tag="cc1type0";
  }
  else if(opt%10==3){
    cut=CCPXcut()+" && fProtonSimPDG == 2212 && fMuonSimPDG == 13 && fXSimPDG==211 && fNeutMode==11 && fEventTypeSim ==1";
    tag="cc11type1";
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  //has to go before TH1D definition
  gStyle->SetOptStat("ou");
  gStyle->SetStatY(0.9);
  gStyle->SetStatX(1);
  gStyle->SetStatW(0.2);
  gStyle->SetStatStyle(0);
  gStyle->SetTitleX(0.5);
  gStyle->SetTitleW(1);

  TH1D * hh = 0x0;

  if(opt<10){
    hh=new TH1D("hh","",2600,-300,2300);
    var="fECalPDGSim";
    tag+="pdg";
  }
  else{
    hh=new TH1D("hh", "", 30, 0,30);
    var="fNeutMode";
    tag+="mode";
  }

  hh->SetXTitle(var);
  
  gPad->SetLogy();

  TString totcut, alltit;

  SetTagTitle(fin, tag, alltit);

  totcut=cut+" && fNECalTracks!=0";
  hh->SetTitle(alltit+totcut); tt->Draw(var+">>hh", totcut);
  cc->Print(Form("outplot/%secal.png", tag.Data()));

  totcut+= "&& fECalPDGSim!=2212 && fECalPDGSim!=2112";
  hh->SetTitle(alltit+totcut); tt->Draw(var+"", totcut);
  cc->Print(Form("outplot/%sNoNucleonecal.png", tag.Data()));

  totcut+= "&& abs(fECalPDGSim)!=321";
  hh->SetTitle(alltit+totcut); tt->Draw(var+"", totcut);
  cc->Print(Form("outplot/%sNoNucleonNoChargeKecal.png", tag.Data()));

  totcut+= "&& abs(fECalPDGSim)!=211";
  hh->SetTitle(alltit+totcut); tt->Draw(var+"", totcut);
  cc->Print(Form("outplot/%sNoNucleonNoChargeKNoChargePiecal.png", tag.Data()));

  totcut+= "&& fECalPDGSim!=130";
  hh->SetTitle(alltit+totcut); tt->Draw(var+"", totcut);
  cc->Print(Form("outplot/%sNoNucleonNoChargeKNoChargePiNoKLongecal.png", tag.Data()));
}

void plotDeltaPTT(const TString fin, const Int_t opt, const Bool_t kmass=kFALSE)
{
  TFile::Open(fin);

  style::SetGlobalStyle();
  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);
  gPad->SetRightMargin(0.05);

  TString tag, cut, truecut, var, alltit;
  if(opt==0){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGModeTypecut(alltit);
    tag="Sim_ModeType";
  }
  else if(opt==1){
    var="fDeltaPTTRec";
    alltit=";#delta p_{TT}^{rec} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGModeTypecut(alltit);
    tag="Rec_ModeType";
  }
  else if(opt==2){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGModeAlltypecut(alltit);
    tag="Sim_Mode";
  }
  else if(opt==3){
    var="fDeltaPTTRec";
    alltit=";#delta p_{TT}^{rec} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGModeAlltypecut(alltit);
    tag="Rec_Mode";
  }
  else if(opt==4){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGBgcut(alltit);
    tag="Sim_BgMode";
  }
  else if(opt==5){
    var="fDeltaPTTRec";
    alltit=";#delta p_{TT}^{rec} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGBgcut(alltit);
    tag="Rec_BgMode";
  }
  else if(opt==6){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCPXcut();
    truecut=CCPXTruePDGModeBreakUpcut(alltit);
    tag="Sim_breakup";
  }
  else if(opt==7){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCMuPPcut();
    truecut=CCMuPPTruePDGTypecut(alltit);
    tag="Sim_MuPP";
  }
  else if(opt==8){
    var="fDeltaPTTSim";
    alltit=";#delta p_{TT}^{sim} (GeV/c)";
    cut=CCPXWrongChargeBgcut(alltit);
    truecut="";
    tag="Sim_WrongCharge";
  }
  else if(opt==20){
    var="fProtonSimDeltaPt->Phi()";
    alltit=";#delta#phi_{T}^{sim} (GeV/c)";
    cut=CCQEcut();
    truecut=CCQETruePDGModeTypecut(alltit);
    tag="CCQESim_ModeType";
  }
  else if(opt==21){
    var="fProtonSimDeltaPt->Phi()";
    alltit=";#delta#phi_{T}^{sim} (GeV/c)";
    cut=CCQEcut();
    truecut=CCQETruePDGModeBreakUpcut(alltit);
    tag="CCQESim_ModeBreakup";
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }
  
  SetTagTitle(fin, tag, alltit);

  TH1D * hh=0x0;
  if(kmass){
    var.ReplaceAll("PTT","2PlusMass");
    alltit.ReplaceAll("#delta p_{TT}","M_{#Delta^{++}}");
    tag+="Mass";

    hh=new TH1D("hh","",200,1,2);
  }
  else if(opt>=20){
    hh=new TH1D("hh","",100,1e-2, 3.142);
    NeutrinoTools::BinLog(hh->GetXaxis());
    gPad->SetLogx();
    gPad->SetLogy();
  }
  else{
   hh=new TH1D("hh","",200,-0.7,0.7);
 }

  hh->SetTitle(Form("%s;p.d.f.", alltit.Data()));
  style::ResetStyle(hh);

  const TString totcut = cut+ truecut;
  const Int_t ent = tt->Draw(var+">>hh",totcut);
  printf("plotDeltaPTT tag %s ent %d var %s totcut { %s }\n", tag.Data(),  ent, var.Data(), totcut.Data());

  //hh->Scale(1./hh->Integral(0,10000)/hh->GetBinWidth(1));
  TH1D * hpdf = NeutrinoTools::ToPDF(hh);
  delete hh;
  hh=hpdf;
  if(gPad->GetLogy()){
    hh->SetMinimum(5e-3);
  }
  else{
    hh->SetMinimum(0);
  }
  hh->SetMaximum(opt%2==0?10:10);
  hh->Draw();

  TLine * ln = new TLine(0,0,0,hh->GetMaximum());
  ln->SetLineColor(kGray);
  ln->SetLineWidth(1);
  ln->Draw("same");
  cc->Print(Form("outplot/%s.png", tag.Data()));
}

void plot2DkT(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();
  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);
  gPad->SetRightMargin(0.05);

  TString tag, cut, truecut, var,alltit;
  if(opt==0){
    var="fProtonSimKt:fXSimKt";
    cut=CCPXcut();
    alltit=";p_{TT}^{#pi+, sim} (GeV/c);p_{TT}^{p, sim} (GeV/c)";
    truecut=CCPXTruePDGModeAlltypecut(alltit);
    tag="Sim2D_Mode";
  }
  else if(opt==1){
    var="fProtonSimKt:fXSimKt";
    cut=CCPXcut();
    alltit=";p_{TT}^{#pi+, sim} (GeV/c);p_{TT}^{p, sim} (GeV/c)";
    truecut=CCPXTruePDGBgcut(alltit);
    tag="Sim2D_bgMode";
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  SetTagTitle(fin, tag, alltit);
  
  TH2D * hh=new TH2D("hh","",500,-0.5,0.5,500,-0.5,0.5);
  style::ResetStyle(hh);
  hh->SetTitle(alltit);
  const Int_t ent = tt->Draw(var+">>hh",cut+ truecut);
  printf("plot2DkT tag %s ent %d\n", tag.Data(),  ent);

  cc->Print(Form("outplot/%s.png", tag.Data()));

}

void plotCCQEvsMuonPt(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();
  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);
  gPad->SetRightMargin(0.15);

  TString tag, cut, truecut, var,alltit;
  Double_t ymin=-999, ymax=-999;
  Double_t hmin=-999, hmax = -999;
  Bool_t klogy;
  Int_t nybin = 20;
  if(opt%10==0){
    var="fProtonSimDeltaPt->Mag():fMuonSimPt->Mag()";
    cut=CCQEcut();
    alltit=";p_{T}^{#mu, sim} (GeV/c);#deltap_{T}^{sim} (GeV/c)";
    tag="ccqedptvsMuonPt";
    ymin = 0;
    ymax = 1;
    hmin = 0;
    hmax = opt<10?0.3:0.6;
    klogy=kFALSE;
  }
  else if(opt%10==1){
    var="fProtonSimDeltaPt->Phi():fMuonSimPt->Mag()";
    cut=CCQEcut();
    alltit=";p_{T}^{#mu, sim} (GeV/c);#delta#phi_{T}^{sim} (rad.)";
    tag="ccqedPHItvsMuonPt";
    ymin = 1e-2;//0;
    ymax = 3.142;
    hmin = ymin;
    hmax = 1;//ymax;//
    klogy=kTRUE;
  }
  else if(opt%10==2){
    var="fProtonSimDeltaPt->Theta():fMuonSimPt->Mag()";
    cut=CCQEcut();
    alltit=";p_{T}^{#mu, sim} (GeV/c);#delta#alpha_{T}^{sim} (rad.)";
    tag="ccqedALPHAtvsMuonPt";
    ymin = 0;
    ymax = 3.142;
    hmin = ymin;
    hmax = 3.142;
    klogy=kFALSE;
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  if(opt<10){
    truecut=CCQETruePDGModeTypecut(alltit);
    tag+="NoEm";
  }
  else{
    truecut=CCQETruePDGModeBreakUpcut(alltit);
    tag+="WithEm";
  }

  SetTagTitle(fin, tag, alltit);
  
  TList * ll= new TList;

  TH2D * hh=new TH2D("hh","",15, 0.1, 2, nybin, ymin, ymax); ll->Add(hh);
  NeutrinoTools::BinLog(hh->GetXaxis());
  if(klogy){
    NeutrinoTools::BinLog(hh->GetYaxis());
  }
  style::ResetStyle(hh);
  hh->SetTitle(alltit);
  const TString totcut = cut+ truecut;
  const Int_t ent = tt->Draw(var+">>hh", totcut);
  printf("plot2DkT tag %s totcut { %s } ent %d\n", tag.Data(),  totcut.Data(), ent);

  TH1D * hnor = 0x0;
  TH1D * hmpv = 0x0;
  TH1D * hwid = 0x0;
  TH1D * hres = 0x0;
  TH1D * hchi = 0x0;
 
  NeutrinoTools::FitSlicesY(hh, hnor, hmpv, hwid, hres, hchi, "RMS", 20, ll);

  //hh->RebinY(1600);

  gPad->SetLogx();
  gPad->SetLogy(klogy);
  hh->Draw("colz");
  cc->Print(Form("outplot/%s.png", tag.Data()));

  TH2D * h2 = NeutrinoTools::NormalHist(hh, 20, kTRUE); ll->Add(h2);
  h2->Draw("colz");
  cc->Print(Form("outplot/%s2dnor.png", tag.Data()));

  //=====================
  gPad->SetRightMargin(0.05);
  //gPad->SetLogy(klogy);
  gPad->SetLogy(0);

  hmpv->SetMinimum(hmin);
  hmpv->SetMaximum(hmax);
  hmpv->Draw();
  cc->Print(Form("outplot/%sMpv.png", tag.Data()));

  hwid->SetMinimum(hmin);
  hwid->SetMaximum(hmax);
  hwid->Draw();
  cc->Print(Form("outplot/%sWid.png", tag.Data()));

  TFile * fout = new TFile(Form("outplot/%sOutFile.root", tag.Data()), "recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

void plotRatioNuclearEmission(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();
  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);
  gPad->SetRightMargin(0.05);

  TString tag, cut, truecut, var,alltit;
  Bool_t klogx;
  Int_t nbin;
  Double_t xmin, xmax, gmax=10, lx;
  if(opt==0){
    var="fProtonSimDeltaPt->Phi()";
    cut=CCQEcut();
    alltit=";#delta#phi_{T}^{sim} (rad.);p.d.f.";
    truecut=CCQETruePDGModeAlltypecut(alltit);
    tag="NuclearEmissionPHI";
    klogx=1;
    nbin = 20;
    xmin = 1e-2;
    xmax = 3.142;
    gmax = 10;
    lx=0.6;
  }
  else if(opt==1){
    var="fProtonSimDeltaPt->Theta()";
    cut=CCQEcut();
    alltit=";#delta#alpha_{T}^{sim} (rad.);p.d.f.";
    truecut=CCQETruePDGModeAlltypecut(alltit);
    tag="NuclearEmissionALPHA";
    klogx=0;
    nbin = 20;
    xmin = 0;
    xmax = 3.142;
    gmax = 1.05;
    lx=0.3;
  }
  else if(opt==2){
    var="fProtonSimDeltaPt->Mag()";
    cut=CCQEcut();
    alltit=";#deltap_{T}^{sim} (GeV/c);p.d.f.";
    truecut=CCQETruePDGModeAlltypecut(alltit);
    tag="NuclearEmissiondpt";
    klogx=0;
    nbin = 20;
    xmin = 0;
    xmax = 0.7;
    gmax = 7;
    lx=0.6;
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  SetTagTitle(fin, tag, alltit);
  
  //====
  TH1D *h0 = new TH1D("h0","",nbin, xmin, xmax);
  if(klogx) NeutrinoTools::BinLog(h0->GetXaxis());
  style::ResetStyle(h0);
  h0->SetTitle(alltit);

  const TString cut0=cut+ truecut + "&& fEventTypeSim==0";
  const Int_t ent0 = tt->Draw(var+">>h0",cut0);
  printf("plot2DkT tag %s cut0 { %s } ent %d\n", tag.Data(), cut0.Data(),  ent0);

  //====

  TH1D *h1 = new TH1D("h1","",nbin, xmin, xmax);
  if(klogx) NeutrinoTools::BinLog(h1->GetXaxis());
  style::ResetStyle(h1);
  h1->SetTitle(alltit);

  const TString cut1=cut+ truecut + "&& fEventTypeSim!=0";
  const Int_t ent1 = tt->Draw(var+">>h1",cut1);
  printf("plot2DkT tag %s cut1 { %s } ent %d\n", tag.Data(), cut1.Data(),  ent1);

  //====

  TH1D * g0 = NeutrinoTools::ToPDF(h0);
  TH1D * g1 = NeutrinoTools::ToPDF(h1);

  TLegend * lg = new TLegend(lx,0.75,lx+0.32,0.88);
  style::ResetStyle(lg);
  lg->SetBorderSize(1);

  gPad->SetLogx(klogx);

  g0->SetMinimum(0);
  g0->SetMaximum(gmax);
  g0->SetLineColor(kBlack);
  //g0->SetLineStyle(kDashDotted);
  g0->SetMarkerStyle(20);
  g0->SetMarkerSize(1);
  g0->SetMarkerColor(kBlack);
  g0->Draw("hist e");
  lg->AddEntry(g0,"no nuclear emission","lp");

  g1->SetLineColor(kRed);
  //g1->SetLineStyle(kSolid);
  g1->SetMarkerStyle(22);
  g1->SetMarkerSize(1);
  g1->SetMarkerColor(kRed);
  g1->Draw("hist e same");
  lg->AddEntry(g1,"w/ nuclear emission","lp");

  lg->Draw();

  cc->Print(Form("outplot/%s.png", tag.Data()));

  NeutrinoTools::ScaleToRef(g1,g0);
  NeutrinoTools::ScaleToRef(g0,g0);
  g1->SetMinimum(0);
  g1->SetMaximum(2);
  g1->Draw("hist e");
  g0->Draw("hist e same");
  lg->Draw();
  g1->GetYaxis()->SetTitle("scaled to ref.");
  cc->Print(Form("outplot/%sratio.png", tag.Data()));

}


void plotEmissionFraction(const TString fin, const Int_t opt)
{
  TFile::Open(fin);

  style::SetGlobalStyle();
  gStyle->SetOptTitle(1);

  TTree *tt = (TTree*)gDirectory->Get("tree");

  TCanvas *cc=new TCanvas;
  style::PadSetup(cc);
  gPad->SetRightMargin(0.05);

  TString tag, cut, truecut, var,alltit;
  Bool_t klogx;
  Int_t nbin;
  Double_t xmin, xmax, gmax=10, lx;
  if(opt==0){
    var="fMuonSimPt->Mag()";
    cut=CCQEcut();
    alltit=";p_{T}^{#mu,sim} (GeV/c);p.d.f.";
    truecut=CCQETruePDGModeAlltypecut(alltit);
    tag="emissionFraction";
    klogx=1;
    nbin = 20;
    xmin = 0.1;
    xmax = 2;
    gmax = 3.5;
    lx=0.6;
  }
  else if(opt==1){
    var="fMuonSimPt->Mag()";
    cut=CCQEcut();
    alltit=";p_{T}^{#mu,sim} (GeV/c);p.d.f.";
    truecut=CCQETruePDGModeAlltypeNoECalcut(alltit);
    tag="NOECalemissionFraction";
    klogx=1;
    nbin = 20;
    xmin = 0.1;
    xmax = 2;
    gmax = 3.5;
    lx=0.6;
  }
  else{
    printf("wrong opt %d\n", opt); exit(1);
  }

  SetTagTitle(fin, tag, alltit);
  
  //====
  TH1D *h0 = new TH1D("h0","",nbin, xmin, xmax);
  if(klogx) NeutrinoTools::BinLog(h0->GetXaxis());
  style::ResetStyle(h0);
  h0->SetTitle(alltit);

  const TString cut0=cut+ truecut;
  const Int_t ent0 = tt->Draw(var+">>h0",cut0);
  printf("plot2DkT tag %s cut0 { %s } ent %d\n", tag.Data(), cut0.Data(),  ent0);

  //====

  TH1D *h1 = new TH1D("h1","",nbin, xmin, xmax);
  if(klogx) NeutrinoTools::BinLog(h1->GetXaxis());
  style::ResetStyle(h1);
  h1->SetTitle(alltit);

  const TString cut1=cut+ truecut + "&& fEventTypeSim==0";
  const Int_t ent1 = tt->Draw(var+">>h1",cut1);
  printf("plot2DkT tag %s cut1 { %s } ent %d\n", tag.Data(), cut1.Data(),  ent1);

  //====

  TH1D * g0 = NeutrinoTools::ToPDF(h0);
  TH1D * g1 = NeutrinoTools::ToPDF(h1);

  TLegend * lg = new TLegend(lx,0.75,lx+0.32,0.88);
  style::ResetStyle(lg);
  lg->SetBorderSize(1);

  gPad->SetLogx(klogx);

  g0->SetMinimum(0);
  g0->SetMaximum(gmax);
  g0->SetLineColor(kBlack);
  //g0->SetLineStyle(kDashDotted);
  g0->SetMarkerStyle(20);
  g0->SetMarkerSize(1);
  g0->SetMarkerColor(kBlack);
  g0->Draw("hist e");
  lg->AddEntry(g0,"all","lp");

  g1->SetLineColor(kRed);
  //g1->SetLineStyle(kSolid);
  g1->SetMarkerStyle(22);
  g1->SetMarkerSize(1);
  g1->SetMarkerColor(kRed);
  g1->Draw("hist e same");
  lg->AddEntry(g1,"no nuclear emission","lp");

  lg->Draw();

  cc->Print(Form("outplot/%s.png", tag.Data()));

  h0->Sumw2();
  h1->Sumw2();
  h1->Divide(h1,h0,1,1,"B");
  h1->SetMinimum(0);
  h1->SetMaximum(1.1);
  h1->Draw("hist e");
  h1->GetYaxis()->SetTitle("fraction of no nuclear emission");
  cc->Print(Form("outplot/%sFraction.png", tag.Data()));

}

